
/****************************************************************************/
/*Copyright (c) 2011, Florent DEVILLE.                                      */
/*All rights reserved.                                                      */
/*                                                                          */
/*Redistribution and use in source and binary forms, with or without        */
/*modification, are permitted provided that the following conditions        */
/*are met:                                                                  */
/*                                                                          */
/* - Redistributions of source code must retain the above copyright         */
/*notice, this list of conditions and the following disclaimer.             */
/* - Redistributions in binary form must reproduce the above                */
/*copyright notice, this list of conditions and the following               */
/*disclaimer in the documentation and/or other materials provided           */
/*with the distribution.                                                    */
/* - The names of its contributors cannot be used to endorse or promote     */
/*products derived from this software without specific prior written        */
/*permission.                                                               */
/* - The source code cannot be used for commercial purposes without         */ 
/*its contributors' permission.                                             */
/*                                                                          */
/*THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS       */
/*"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT         */
/*LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS         */
/*FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE            */
/*COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,       */
/*INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,      */
/*BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;          */
/*LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER          */
/*CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT        */
/*LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN         */
/*ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE           */
/*POSSIBILITY OF SUCH DAMAGE.                                               */
/****************************************************************************/

#include "RTPrimTore.h"

#include <math.h>
#include "RTPolynomial.h"

namespace RT
{
	/*Constructor*/
	RTPrimTore::RTPrimTore(F32 bigRadius, F32 smallRadius):RTIPrimitive(), m_bigRadius(bigRadius), 
			m_smallRadius(smallRadius){}

	/*Destructor*/
	RTPrimTore::~RTPrimTore(){}

	/*Calculate the intersection between a ray an the torus*/
	bool RTPrimTore::intersect(const RTRay& ray, F32& dist, bool useGlobalCoordinates)
	{
		//change the coordinate system
		RTRay R;
		if (useGlobalCoordinates)
			globalToLocal(ray, R);
		else
			convertRayToLocalCoordinates(ray, R);

		F32 alpha = R.getDirection().dot(R.getDirection());
		F32 beta = 2 * R.getDirection().dot(R.getOrigin());
		F32 gamma = R.getOrigin().dot(R.getOrigin()) - m_smallRadius * m_smallRadius - m_bigRadius * m_bigRadius;

		F64 coef[5];
		coef[4] = alpha * alpha;
		coef[3] = 2 * alpha * beta;
		coef[2] = beta * beta + 2 * alpha * gamma + 4 * m_bigRadius * m_bigRadius * R.getDirection().z * 
			R.getDirection().z;
		coef[1] = 2 * beta * gamma + 8 * m_bigRadius * m_bigRadius * R.getOrigin().z * R.getDirection().z;
		coef[0] = gamma * gamma + 4 * m_bigRadius * m_bigRadius * R.getOrigin().z * R.getOrigin().z - 
			4 * m_bigRadius * m_bigRadius * m_smallRadius * m_smallRadius;

		//create the RTPolynomial
		RTPolynomial p(coef[4], coef[3], coef[2], coef[1], coef[0]);

		//check the degree and gets it
		p.Check();
		//I32 n = p.Degree();

		//solve the equation
		F64 solutions[4] = {0, 0, 0, 0};
		p.SolveQuartic(solutions);

		//get the greatest solution
		F64 min = -1;//solutions[0];
		for(I32 i = 0; i < 4; i++)
		{
			if(solutions[i] > min)
				min = solutions[i];
		}

		if(min <= 0)
			return false;

		//get the smallest solution
		for(I32 i = 0; i < 4; i++)
		{
			if(solutions[i] < min && solutions[i] > 0)
				min = solutions[i];
		}

		if (min <= 0)
			return false;

		//set the distance
		dist = (F32)min;
		return true;

	}

	/* 
	To calculate the normal, we use the gradient
	Equation : (x^2 + y^2 + z^2 + R^2 - r^2)^2 - 4R^2(x^2 + z^2)
	Gradiant : 
				df/dx = 4x(x^2 + y^2 + z^2 + R^2 - r^2) - 8xR^2
				df/dy = 4y(x^2 + y^2 + z^2 + R^2 - r^2)
				df/dz = 4z(x^2 + y^2 + z^2 + R^2 - r^2) - 8zR^2
	Simplify by 4:
				df/dx = x(x^2 + y^2 + z^2 + R^2 - r^2) - 2xR^2
				df/dy = y(x^2 + y^2 + z^2 + R^2 - r^2)
				df/dz = z(x^2 + y^2 + z^2 + R^2 - r^2) - 2zR^2
	*/
	void RTPrimTore::computeNormal(const RTRay& ray, F32 dist, RTVector3f& normal)const
	{
		//change the coordinate system
		RTRay R;
		convertRayToLocalCoordinates(ray, R);

		//calculate intersection poI32
		RTVector3f i = R.getOrigin() + (R.getDirection() * dist);
	
		//temp variable to calculate the gradient
		F32 R2 = m_bigRadius * m_bigRadius;
		F32 r2 = m_smallRadius * m_smallRadius;
		F32 n = i.norme2() - R2 - r2;

		//use the gradient
		normal.x = 4 * i.x * n;
		normal.y = 4 * i.y * n; 
		normal.z = (4 * i.z * n) + (8 * R2 * i.z);
		normal.normalize();

		//if inside
		F32 d = R.getDirection().dot(normal);
		if(d > 0)
			normal = normal * -1;

		//reset to the global coordinate systeme
		convertNormalToGlobalCoordinates(normal, normal);	
	}

	/*Check if the point is inside the torus*/
	bool RTPrimTore::isInside(const RTVector3f& p)const
	{
		//convert p to local coordinates
		RTVector3f localP;
		convertPointToLocalCoordinates(p, localP);

		F64 p1 = localP.norme2() + m_bigRadius * m_bigRadius - m_smallRadius * m_smallRadius;
		F64 p2 = 4 * m_bigRadius * m_bigRadius * (localP.x * localP.x + localP.z * localP.z);

		if(p1 <= p2)
			return true;

		return false;
	}

}
